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Abstract -Isostatic networks are minimally rigid and therefore have, generically, nonzero elastic 
moduli. Regular isostatic networks have finite moduli in the limit of large sizes. However, numer- 
ical simulations show that all elastic moduli of geometrically disordered isostatic networks go to 
zero with system size. This holds true for positional as well as for topological disorder. In most 
cases, elastic moduli decrease as inverse power-laws of system size. On directed isostatic networks, 
however, of which the square and cubic lattices are particular cases, the decrease of the moduli is 
exponential with size. For these, the observed elastic weakening can be quantitatively described 
in terms of the multiplicative growth of stresses with system size, giving rise to bulk and shear 
moduli of order e~ bL . The case of sphere packings, which only accept compressive contact forces, 
is considered separately. It is argued that these have a finite bulk modulus because of specific 
correlations in contact disorder, introduced by the constraint of compressivity. We discuss why 
their shear modulus, nevertheless, is again zero for large sizes. A quantitative model is proposed 
that describes the numerically measured shear modulus, both as a function of the loading angle 
and system size. 

In all cases, if a density p > of overconstraints is present, as when a packing is deformed by 
compression or when a glass is outside its isostatic composition window, all asymptotic moduli 
become finite. For square networks with periodic boundary conditions, these are of order y/p. 
For directed networks, elastic moduli are of order e _C//p , indicating the existence of an "isostatic 
length scale" of order 1/p. 
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Disorder in condensed matter gives rise to interesting 
physics that has been the object of intense study for 
decades [TJ Important problems to which much ef- 
fort is still being devoted are the vibrations and elec- 
tronic conduction in disordered solids pQ. Static prop- 
erties like the elastic moduli of disordered systems, on the 
other hand, appear to be more readily accessible to anal- 
ysis, so one would expect them to hold fewer surprises. 
This expectation turns out to be true only for hyperstatic 
networks. Consider, for example, a triangular spring net- 
work whose sites have been randomly displaced within a 
circle of radius e, with repose lengths of all springs ad- 
justed accordingly. All elastic moduli of these positionally 
disordered (PD) networks are finite and size-independent 
for large sizes, and go continuously to their homogeneous 
limit when e — > [3] . Formally, Yl (e) sa — ae 2 for large 
L. A similar behavior is observed for topologically disor- 
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dered (TD) networks obtained e.g. by locally rewiring a 
small fraction e of the links. 

However, the simple perturbative picture described above 
is not universal. It is only valid for elastic networks that, 
like the triangular network, are Statically Indeterminate or 
Hyperstatic, i.e. have more springs than strictly required 
to be rigid. Networks which are minimally rigid and do 
not have redundant springs are Statically Determinate or 
Isostatic. Their clastic properties show, as their size in- 
creases, a dramatic disorder-weakening effect that is the 
subject of this work. We will argue that all elastic moduli 
of isostatic networks go to zero with system size, if these 
have random geometrical (positional or topological) disor- 
der. 

An Isostatic Network (IN) is mechanically rigid, and 
therefore has no flexibilities (zero-energy modes), other 
than the trivial m(d) = d(d + l)/2 rigid roto-translations 
in d dimensions. But the removal of any spring introduces 
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Fig. 1: a) Random isostatic network with periodic boundaries. 
These networks are built by means of a matching algorithm for 
rigidity [4]. b) Positionally disordered square lattice attached 
to a fixed boundary (black circles) and loaded at the top, with 
periodic boundary condition in the horizontal direction. 



one flexibility, that is, an isostatic network has no redun- 
dant springs. A simple counting of degrees of freedom (d 
degrees of freedom per site) and constraints (each spring 
counts as a constraint), originally due to Maxwell [5], 
shows that a freely-standing isostatic network of N points 
in d dimensions has exactly dN — m(d) springs. This 
condition is called constraint balance. Neglecting bound- 
ary contributions, their average coordination is therefore 
< z >— 2d. It is important to remark that constraint 
balance is not a sufficient condition for isostaticity. Addi- 
tionally, all springs must be "properly distributed" , in the 
sense defined in the theory of graph rigidity. 

While building topologically disordered isostatic net- 
works is in general hard (satisfying constraint balance is 
easy, verifying proper distribution is highly nontrivial), 
in two dimensions simple integer algorithms exist [4 . 6 
for this task. On the other hand, topologically regular 
isostatic lattices are easily obtained. The simplest 
examples of such regular INs are the square and cubic 
lattices with appropriate boundary conditions. These 
are topologically ordered, as every site in the bulk has 
coordination z = 2d. Given that square and cubic lattices 
are so common, it is surprising that their elastic proper- 
ties in the presence of positional disorder have not been 
discussed in detail before. Furthermore, the relevance of 
our findings is not limited to topologically regular iso- 
static networks. In recent years, it has been noticed that 
several natural systems self-organize onto topologically 
disordered structures that are isostatic, i.e. minimally 
rigid without redundant links. Following the discovery of 
isostaticity [7] in frictionless sphere-packings (which are 
models of metallic glasses [8], colloids [9], emulsions [TP] , 
and granular matter [Tllllj). it was later proposed [12, 13 
that covalent glasses are isostatic in a certain composition 
window. These further instances of isostatic networks 
are topologically disordered, as the coordination number 
fluctuates from site to site. Isostatic systems have a 
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Fig. 2: a) Bulk elastic modulus Bl{p) for PDSQ with pe- 
riodic boundary conditions (PBC) in the isostatic case (full 
circles) and with a density p of excess springs equal to 1 x 10 -4 
(plusses), 5xl0~ 4 (crosses), 2xl0 -3 (stars), lxlCT 2 (squares), 
and 5 x 10~ 2 (circles). Solid lines are fits to the data (see text), 
b) Asymptotic modulus B ao (p) vs p for PDSQ with PBC. The 
full line is a power law fit with exponent 7 = 0.53 ± 0.05. 



number of anomalous properties that derive from their 
being at the verge of loosing mechanical stability. Much 
attention has been devoted in recent years to the study 
of the particular static [71 IT4HT9] and vibrational [T8ti20] 
properties that isostatic systems display. 
In this work, elastic properties arc explored for two 
different types of disordered isostatic networks. These 
are: positionally disordered square networks (PDSQ), 
and topologically disordered isostatic networks (TD). 
The latter are built as a random subset of the triangular 
lattice, enforcing isostaticity by means of matching 
algorithm for rigidity [1]. An example of TD is shown 
in Fig. QJ,. For the purpose of these studies, e = 0.2 is 
taken for all disordered networks. Elastic properties of 
isostatic networks are strongly dependent on boundary 
conditions. Let us first consider two-dimensional isostatic 
networks with periodic boundary conditions, made of 
N = L 2 sites with 2N (properly distributed) springs 0. 



1 Since two translations are isometrics on a torus, two of these 
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Fig. 3: a) Bulk elastic modulus Bl(p) for TD networks with 
PBC in the isostatic case (full circles) and with a density p of 
excess springs equal to 1 x 10~ 3 (plusses), 5 x 10 -3 (crosses), 
2 x 10~ 2 (stars), 1 x 10" 1 (squares), and 5 x 10" 1 (circles). 
Solid lines are fits (see text), b) Asymptotic modulus B oa (p) 
vs p for TD networks with PBC. 



The bulk elastic modulus B of an elastic network with 
periodic boundary conditions can be measured as the 
per-link elastic energy due to an infinitesimal expansion 
of all springs. This modulus is finite and size independent 
for regular (topologically as well as positionally ordered) 
networks as e.g a regular square lattice. For PDSQ and 
TD networks, on the other hand, there is a disorder- 
induced elastic collapse, whereby B 1 ^ goes to zero as a 
power law of system size L, as shown in Figs. [2jt and [3^,. 
From fits to these data it is found that -B]f ~ aL~^ with 
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/r = 1.3 ± 0.1, and m TD = 5.5 ± 0.5. The value of 

\l is not universal but dependent on e, network structure, 
and boundary conditions. Furthermore, as we later show, 
in certain cases the decay of the modulus is not even 
power-law but exponential with size. 
This disorder-induced elastic collapse only happens for 



springs are redundant. Therefore these networks are not strictly 
isostatic. However, their number of redundancies is not extensive 
but small, and can be ignored for the purpose of our discussion. We 
will refer to them as being isostatic. 



isostatic networks. If a small density p > of extra 
springs per site is now added [2T] , these networks become 
overconstrained (or hyperstatic) with average coordina- 
tion z = 4 + 2p. Overconstrained disordered networks 
have, in the large-size limit (LSL), a finite bulk modulus 
Boo{p) (See Figs.rjb and[3jD). In order to estimate B^p) 
for asymptotically large sizes, we fit the bulk modulus 
data for p > using ad-hoc expressions, since at present 
there is no theoretical prediction for the functional form 
of finite-size effects. For PDSQ networks with p > 
we use B L {p) = Bf e~ c{p)L + 5oo(»(l - e~ c ^ L ), with 
c(p) and i?oo(p) as fitting parameters for each p. From 
the resulting fits we find that c(p) ~ bp a with a « 0.5. 
This behavior would suggest the existence of a crossover 
length of order p" 1 ! 2 . We notice that a characteristic 
length-scale in this problem is the average distance be- 
tween overconstrained links, and that this length behaves 
precisely as p~ x l 2 in two dimensions. The asymptotic 
modulus for PDSQ is found to behave as B oa (p) ~ dp 1 ', 
with 7 = 0.53 ± 0.05 (Fig. rjb). 

For TD networks, a simple expression 
B L {p) = {{Bf°f + (Boo(;p)) e } 1/e , with e(p) and B^v) 
as fitting parameters, is found to describe the data in 
Fig. [5^ acceptably well, and e(p) turns out to be close 
to one in most cases. Our results for B ao (p) of TD 
networks is shown in Fig. [5Jd, where it is seen that no 
single power-law of p describes the data well. 
We now consider the case of directed INs; those for which 
rigidity stems from a boundary, and whose structure 
is such that all spring forces can be obtained by local 
propagation of the loads using force equilibrium. We find 
that, while disorder- induced collapse happens for directed 
isostatic networks as well, its functional dependence 
with system size is entirely different than in the case of 
periodic boundaries. Two examples of disordered directed 
INs discussed here are the square positionally disordered 
directed network (PDDN - An example is shown in 
Fig- Hp) and the topologically disordered directed net- 
work (TDDN). TDDN are built by letting each site be 
supported by two springs, randomly chosen among its 
three lower neighbors on a triangular lattice (For details, 
see jTUTlHT?]). Both PDDN and TDDN arc attached to 
a rigid boundary at the bottom, have periodic boundary 
conditions in the horizontal direction, and a free upper 
boundary where loads are applied. 
Elastic moduli of these networks are proportional to the 
inverse per-spring elastic energy resulting from letting a 
unit load act on each site of their free upper boundary. 
If these loads are vertical, the uniaxial compression 
modulus Y is obtained, while if the loads are horizontal, 
the shear modulus G is measured. In the absence of 
disorder, directed square networks have finite uniaxial 
compression modulus Y and shear modulus G. For 
PDDN and TDDN, on the other hand, both Y and G 
go to zero exponentially with size. Fig. [4^ shows the 
uniaxial compression modulus of PDDN, while Fig. 2b 
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Fig. 4: Uniaxial compressive elastic modulus Y for: a) PDDN 
with e — 0.2, and b) TDDN, with a density of extra springs 
p = (isostatic case, plusses), p — 1CP 3 (crosses), 2 x 10 -3 
(stars), 5 x 10~ 3 (squares), 10~ 2 (circles), 2 x 10~ 2 (triangles), 
5 x 10~ 2 (downwards triangles), and 10" 1 (diamonds). Lines 
are fits using Eq. {TJ. 

shows data for TDDN. The shear modulus G (not shown) 
for PDDN and TDDN behaves in a similar way, going to 
zero exponentially with size. 

Also in the case of directed networks, disorder-induced 
elastic collapse is a consequence of isostaticity. Adding a 
density p > of extra springs at each site, overconstrained 
networks are obtained. These have, in the LSL, a finite 
modulus Yoo(p), which we estimate by fitting (the use of 
this expression will be justified later) 

to the data in Figs. |4jt,b, with Y ao (p), a, and b as free 
parameters when p > 0. For p — >• 0, one has — > and 
dU reduces to Y(L) = abL/(e bL - 1), with only two fitting 
parameters. 

Let us now discuss these observations for directed net- 
works. For an clastic network under a load a per unit 
area, Y oc a 2 / f2 [22] - where f 2 is the second moment of 
the spring-force distribution. The vanishing of Y with in- 
creasing L then implies f 2 — > oc in that limit. Fig. [SJi 
shows the second moment (d) of spring stresses at 



a fixed depth d (distance from the top) for several val- 
ues of p and for L = 100 and 70, on e = 0.20 PDDNs. 
When p = 0, f 2 grows exponentially with depth as f 2 e bd . 
For overconstrained networks with p > 0, the growth of 
f 2 is initially exponential but saturates at a p-dcpcndcnt 
value f 2 m (p). The exponential growth of f 2 for p = can 
be understood in terms of previously reported [71 ll4f[T7] 
multiplicative properties of stress propagation on directed 
INs. We can now provide a justification for Eq. ([1]). For 
PDDN and TDDN, since they are directed, stresses are de- 
termined propagatively from the top down, without any 
information at all coming from below [71 [T4Tl~51[l~7] . There- 
fore f 2 (d) cannot depend on L. This is verified in 
Fig. [5k. with p = 0, where data for L = 100 and 70 
are seen to behave similarly. When p > 0, propagativity 
no longer holds, as stresses on overconstrained networks 
depend on the rheology of the whole system. Despite 
this, Fig. [Sk shows that (d) remains approximately 

L-indepcndent for small p as well, and therefore can be 
approximately described by an L-indcpcndcnt expression. 
We choose f 2 (d) ~ (e- bd /f° + l/if™)" 1 , where f l 2 lm {p) 
is the asymptotic value of f% at large depths. This simple 
ansatz allows us to calculate y _1 (L,p) ~ 1/L J Q L dtf2(t) 
from which Eq. ([T]) results, with Y^ oc l/fj m . As seen m 
Fig.[§L,b, Eq. © fits our data for PDDN and TDDN quite 
well. The resulting asymptotic values Foo(p) for PDDN 
and TDDN are found to be clearly incompatible with a 
power-law form Y OQ (p) ~ p e , but they can be fitted quite 
closely using Y^ = yoe~ c ^ p , as shown in Fig. [SJd. Writing 
Yqo ~ I//2™ ~ e~ bds , where d s {p) is a stress-saturation 
depth, this last result then implies that d s ~ 1/p. The fol- 
lowing picture then emerges: on directed INs, stress fluctu- 
ations grow multiplicatively [7lTl"4"HT7] , their scale increases 
exponentially with depth, and therefore Yl ~ e~ bL . For 
slightly overconstrained networks, this behavior is essen- 
tially unchanged, i.e. the system behaves as if isostatic, 
until a scale d s ~ 1/p is reached, beyond which f 2 satu- 
rates and Yl attains a finite value Y^ . We notice that an 
'isostaticity length' £ ^ 1/p has been defined in the context 
of vibrational modes of isostatic systems [18] . Our results 
show that this same scale determines the elastic proper- 
ties of directed INs. Given that 1/p also determines the 
average vertical separation among redundant springs, this 
might suggest that directed isostatic networks behave as 
if effectively one-dimensional, probably due to unidimen- 
sional character of stress propagation on these systems. 
Further work is needed to verify this picture in three di- 
mensions. 

We remark that disorder-induced elastic collapse happens 
for isostatic networks with any type of local structure and 
any sort of boundary conditions, if these have uncorre- 
cted geometric disorder. Other forms of randomness, like 
mass or spring-constant disorder, which determine the dy- 
namical properties of networks pQ, do not have an elastic 
weakening effect on INs. Additionally, notice that appro- 
priately chosen correlations in geometrical disorder may 
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Fig. 5: a) Second moment /2 of spring stresses at depth d on 
PDDNs with disorder e = 0.20, for three different values of the 
density p of excess springs. Data shown are for sizes L = 100 
(empty symbols) and 70 (filled symbols). TDDNs show a sim- 
ilar behavior, b) Asymptotic compressive modulus Yoo(p) ver- 
sus overconstraint density p, for PDDN (squares) and TDDN 
(circles). Solid lines are fits of the form log Yoo/j/o = —c/p. 



effectively suppress elastic collapse. This is the case with 
compressive isostatic packings, i.e. those built under the 
constraint that no contact be subject to traction, as we 
discuss next. A system of frictionless stiff spheres self or- 
ganizes to satisfy this constraint by rearranging its con- 
tacts, for a given external load (which we call the prepa- 
ration load), by opening contacts subject to traction, and 
replacing them with compressive ones. The resulting iso- 
static set of contacts is disordered, but not independently 
so. Correlations among them have been imposed by the 
comprcssivity condition, under the given preparation load. 
The absence of negative forces makes large stresses impos- 
sible, as we discuss later. As a result, the elastic modulus, 
measured in the direction of the preparation load, is finite 
in the LSL. But for any load other than the preparation 
load (if further contact rearrangements arc not allowed) 
the elastic response vanishes in the LSL. This has been 
observed in numerical studies of compressive sphere pack- 
ings, where the bulk modulus is found to be finite for lin- 
ear elastic interactions, but the shear modulus vanishes 



asymptotically, at the isostatic point [TT| . 
The vanishing of the shear modulus while the compressive 
modulus is finite can be illustrated, and discussed quan- 
titatively, with the help of a toy model for compressive 
sphere packings, which we call Compressive-Stress Topo- 
logically Disordered Directed Network (CSTDDN). This 
model is defined as follows. For a given (e.g. vertical) 
preparation load, contact forces are calculated by propa- 
gation from the top, layer by layer. At each site, sup- 
porting leg configurations are chosen at random, but only 
among the two possibilities (out of three in 2d) which en- 
sure that all contact forces below this site are compressive 
(See [7J[T3J[T5J[T7] ). This compressivity condition makes 
the configuration of contacts dependent upon the loads 
applied at the top, and also upon the particular set of 
contacts chosen at all intermediate steps, thereby intro- 
ducing correlations. The reason why these correlated INs 
have finite modulus in the direction of the correlating load 
is simple. The sum of all vertical stress components across 
a given cut must equal the total load applied at the top, 
i.e. contact forces satisfy a "conservation equation" . Since 
negative forces are absent by construction, contact forces 
cannot grow arbitrarily large with size as they do on un- 
corrected isostatic networks, where large negative forces 
can be balanced by large positive ones so as to still sat- 
isfy this conservation rule. Therefore the scale of forces 
remains finite in the LSL, and so does the compressive 
modulus Y (Fig. [5]) . But these correlations arc ineffective 
to limit the scale of contact forces if now the load at the 
top differs in any way from the preparation load, with- 
out allowing for further contact rearrangements. Consider 
the case of a load that deviates by an angle from the 
preparation load. A CSTDDN, "prepared" with a verti- 
cal homogeneous load o~ v , has a depth-independent second 
moment of stresses f 2 y = Py/Yo and a finite modulus Y a . 
If an additional shear a x is now applied (without mod- 
ifying the network topology), the second moment of the 
additional stresses will grow exponentially with depth x 
as fx* ~ o~%.ae . By superposition, the resulting modulus 
Y~L(o- 2 x + o- 2 y )/f Q L (fZ* + f?)dx is 



Y = Yn 



1 + e 2 



l + aY Q Q 2 {e bL -l)/bV 



(2) 



where 9 = a x /a y measures deviation of the measure- 
ment load from the direction of the preparation load. As 
seen in Fig. |6l our numerical data for the elastic mod- 
ulus in this model is described by this expression very 
well. Simple analysis of ^ reveals three regimes that are 
clearly recognizable in Fig. [SJ where Y behaves respec- 
tively as i) Y ~ Yq for O 2 < e~ bL /Y , ii) Y ~ e 2 e~ bL for 



-bL 



/Y < 6 2 < 1, and iii) Y 



-bL 



for 6 > 1. 



We have thus shown that, while regular isostatic networks 
have finite elastic moduli, the elastic constants of geomet- 
rically disordered isostatic networks vanish in the limit 
of large system sizes, although in a manner that depends 
on structure and boundary conditions. INs with periodic 
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boundary conditions have a bulk elastic modulus that goes 
to zero as a power-law of size L with a non-universal ex- 
ponent. For the case of directed INs grounded to a rigid 
bar and loaded at the opposite end, all moduli go to zero 
exponentially with size 

The vanishing of clastic modulus with size has been pre- 




Fig. 6: Elastic modulus of CSTDDNs, a toy model for com- 
pressive sphere packings, measured at an angle O from the 
preparation load, for sizes L = 100 (squares), 200 (circles), 300 
(triangles), and 400 (diamonds). Solid lines are fits to 

viously reported [13] for self-organized isostatic networks. 
These backbones, being the isostatic subset of a random 
rigidity percolation backbone [53], are necessarily fractal. 
Wc have shown that even compact isostatic structures, in- 
cluding the square and cubic lattices, have vanishing mod- 
uli in the largc-sizc-limit if they are geometrically disor- 
dered. 

As a closing remark, we note that models of self- 
organization in glasses [T3][2I] propose that certain calco- 
genides may form isostatic structures in order to lower 
their elastic energy by avoiding self-stress. This argument 
should be modified, in the light of the results reported 
here, when a large enough external load F = aft acts 
on the melt. The total elastic energy per unit volume is 
e-seif + ei oa( i where ei oa & ~ a 2 /Y is maximized for an iso- 
static structure, as we have shown. The minimal energy 
structure would then be determined by a load-dependent 
balance between two opposing effects: reducing self-stress 
by making the structure closer to isostatic, and increasing 
the elastic modulus by making the structure more over- 
constrained. 

* * * 

The author thanks M. Wyart, P. Duxbury and 
M. V. Chubynsky for useful discussions. 
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